Ergodicity of Thermostat Family of Nose— Hoover type 
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One-variable thermostats are studied as a generalization of the Nose-Hoover method which is 
aimed to achieve Gibbs' canonical distribution with conserving the time-reversibility. A condition 
for equations of motion for the system with the thermostats is derived in the form of a partial 
differential equation. Solutions of this equation construct a family of thermostats including the 
Nose-Hoover method as the minimal solution. It is shown that the one-variable thermostat coupled 
with the one-dimensional harmonic oscillator loses its ergodicity with large enough relaxation time. 
The present result suggests that multi-variable thermostats are required to assure the ergodicity 
and to work as heatbath. 
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One of the important issues of recent simulation stud- 
ies is achieving the canonical distribution for the system 
at the desired temperature. Traditional molecular dy- 
namics (MD) simulation has been performed on the basis 
of the Hamiltonian form which gives the microcanonical 
distribution. In the microcanonical system, it is diffi- 
cult to control the temperature since all we can do is 
to set up the initial configuration. Therefore, we need 
canonical MD which is defined as a method to achieve 
the canonical distribution for the system. In addition 
to the above, some properties are also desired; (i) Au- 
tonomous dynamics, i.e., the equations of motion should 
be a closed form and the dynamics should be determin- 
istic; (ii) Time-reversibility; (iii) Ergodicity. 

Many methods are proposed to control temperature in 
MD simulations. The first method controlling tempera- 
ture was proposed by Woodcock [1]. While this method 
is very simple, it is non-autonomous since it involves 
artificial velocity-scaling. The autonomous method is 
proposed on the basis of the variational principle with 
constraint which is refered to the Gaussian thermo- 
stat [2]. This thermostat gives the canonical distribu- 
tion for potential energy with conserving kinetic energy. 
Nose proposed the extended system method which gives 
the canonical distribution for the total energy [3]. This 
method was reformulated to a simple form by Hoover, 
and it is now refered to the Nose-Hoover method [4]. 
The Nose-Hoover method achieves the canonical distri- 
bution for a given system by adding one degree of free- 
dom [5]. The Hamiltonian formulations have been also 
proposed [6, 7] . Recently, Hoover et al. showed that the 
deterministic thermostats can be applied for far-from- 
cquilibrium problems [8] and Kusnezov et al. extended 
the Nose-Hoover dynamics to classical Lie algebras [9]. 

While the Nose-Hoover method is convenient to study 
various isothermal systems, it is found that the method 
sometimes loses its ergodicity, and consequently fails to 



achieve the canonical distribution. In order to improve 
the ergodicity of the Nose-Hoover method, some ex- 
tended methods are proposed [10-12]. In Ref. [12], Kus- 
nezov et al. proposed the general formulation of the ex- 
tended Nose-Hoover method and concluded that two ad- 
ditional degrees of freedom are enough to make a system 
ergodic. However, we have not had the reason why the 
multi-variable thermostat achieves the ergodicity while 
the single- variable ones lose [13]. Therefore, we study 
the ergodicity of general single-variable thermostats in 
order to investigate when and why the system loses its 
ergodicity. In the present Letter, we first derive the con- 
dition which the equations of motion should satisfy to 
achieve Gibbs' canonical distribution and we give the 
general expressions for the one-variable thermostats ex- 
tended from the Nose-Hoover method. Then we show 
that the one-variable thermostat coupled with the one- 
dimensional harmonic oscillator loses its ergodicity for 
large relaxation time. 

Consider the distribution function / and the state vec- 
tor r of a phase space. Let H(T) be a pseudo Hamil- 
tonian describing the energy of the system at T. The 
distribution function is normalized as 



far = i, 



(i) 



and the internal energy of the system is given by U as 

HfdT = U. (2) 
The entropy of the system is defined by 
S = -k B [ /log/dT 



with the Boltzmann constant fee- The equilibrium state 
is obtained by maximizing the entropy under the condi- 
tions Eqs. (1) and (2), and the canonical distribution is 
obtained to be 
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(3) 
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with the partition function Z = J exp (— pH)dT. In or- 
der to perform an MD simulation, equations of motion for 
r must be explicitly given. The equation of continuum 
for / and T is 

^(17) = o, 

where d/dT denotes the divergence and the distribution 
is assumed to be stationary. In order to achieve the 
canonical distribution Eq. (3), we have the following con- 
dition for r as 



0T 



(4) 



Note that the flow of this dynamics is compressible since 
the divergence of T is not zero, while the flow is incom- 
pressible in the microcanonical system [12]. 

The equations of motion satisfying Eq. (4) achieve the 
canonical distribution for arbitrary chosen H, provided 
that the system is ergodic. In most cases, the system of 
interest is described by a Hamiltonian. Let H a be such 
Hamiltonian defined in a 27V-dimensional phase space 
To = (<7ij ■ ■ ■ j (In, Pi, ■ ■ ■ ,Pn) which is a subspace of T, 
that is, T = T ® T±. The distribution function of the 
subsystem is obtained by the projection from T onto To 
as 

/o(T ) = //dT±. (5) 

If the pseudo Hamiltonian H is chosen as 
H(T) = H (T )+H ± (T ± ), 
then the distribution function becomes 



f(T) = f (T )f±(T ± ), 



(6) 



since / oc exp(—/3H). With Eqs. (5) and (6), we obtain 
the canonical distribution for the given Hamiltonian to 
be 

fo = Zo 1 exp(-0H o ) 

with Zq 1 = Z- 1 j exp (-0H±)dT±. 

Even if the pseudo Hamiltonian H is explicitly given, 
there are various choices of the dynamics. In the present 
Letter, we consider one- variable thermostats as extension 
of the Nose-Hoover method since it is favorable to simu- 
late systems with less degrees of freedom. Then the total 
phase space is defined by T = (qi, • • • , qN,Pi, ■ ■ ■ ,Pn) <S> 
(£) with the additional degree of freedom £. In order that 
the distribution function of the subsystem exists, the in- 
tegration of the total distribution function over £ should 
converge as 
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J — < 



exp (-f3H±)d( < oo. 
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FIG. 1: The range of the energy Ho- The function <j>(Ho) = 
Ho — (rn + l)/3~ log Ho and C are plotted for m = 0, (3 = 1, 
and C — 1.019. The minimum and the maximum values 
are determined as solutions of 4>{Ho) = C. Note that this 
equation always has two positive solutions ff m i n and ff mal . 
From the inequality (15), we can estimate the range of the 
energy as 0.815 < H < 1.211. 



The simplest function satisfying the above condition is 
H±(() = ( 2 /2, and therefore we consider the pseudo 
Hamiltonian 



H 



1 



2 /-2 



Ho + ~ 2 rX 



(7) 



with the relaxation time r. For the simplicity, we con- 
sider the subsystem Hq with one degree of freedom here- 
after. The following arguments are not changed in the 
case with many degrees of freedom. 

Consider the following equations of motion: 



dH 
dq 
dH 
dp 



(8) 

(9) 
(10) 



which are simple extension of the Nose-Hoover method. 
The function g(p, Q is a friction term which is p£ in the 
Nose-Hoover method. The time derivative of ( depends 
on g and it is determined from the condition Eq. (4). 
From Eqs. (4) and (7), we have the following partial dif- 
ferential equation: 



which C should satisfy. Here we assumed the natural 
Hamiltonian form Hq = p 2 /2 + V(q) with the potential 
energy V. The function g depends only on p and ( since 
Eq. (11) does not contain q. The solution of the equation 
gives £ as a function of p and £, and then the equations of 
motion are closed and become autonomous. The solution 
of Eq. (11) for the case d(/d( = is given in Ref. [12]. 
In the present Letter, we study more general solutions 
both for d(/d( ^ 0. 
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The equations of motion of the Nose-Hoover method 
are time-reversible with the operation p — p, q — > q, 
and £ — > — C- Similarly, the equations of motion Eqs. (8)- 
(10) are also time- reversible when g — > g. Therefore, the 
function g is a linear combination of p k ( l (k > 0,1 > 
0, k + I = 2, 4, 6, • • •). Here we assume that g does not 
contain the negative power of p and £ for the stability of 
the dynamics. In the case that both k and I are even, 
it is difficult to control temperature since p k ( l becomes 
positive scmi-dchnite [14]. Therefore, we consider only 
odd cases as g = p2m+i^2n+i ( m n = 0, 1,2, •••), and 
then we obtain the expression for £ as 



1.4 



C = ~^Zn P 



2m+2 



2m + 1 







-P 



,2 m 



(12) 



where the function z n (Q is the solution of the following 
ordinary differential equation: 



and it is explicitly expressed as 



2-~i h\\ 9 



P n -^k\ \ 2 



Equation (12) gives the Nose-Hoover method as the 
minimal solution with (to, n) = (0, 0). The case (to, n) = 
(1,0) gives the thermostat controlling only the second 
moment of the kinetic energy {K 2 ) as 



P% 
1 

72 
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The kinetic-moments method [11] is obtained from g = 
9i(p,0 + 92(P,0 with for (m, n) = (0,0) and g 2 for 
(m,n) = (l,0). 

In order to study the ergodicity of the present extended 
method, we consider the one-dimensional harmonic os- 
cillator described by the Hamiltonian H = p 2 /2 + q 2 /2. 
Then we have the following equations of motion: 
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FIG. 2: (Color online) Time evolutions of the energies of the 
systems with three thermostats (a) g — p£, (b) g = p£ , and 
(c) g = p 3 C- Three cases t = 5, 10, and 50 are plotted in each 
figure. The upper and the lower limits of energies estimated 
from the inequality (15) are shown as the solid lines. The 
theoretical estimation becomes more accurate with larger r. 
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Introducing the polar coordinates by p — rcosO and 
q = r sin 9, we have the equations of motion in terms 
of (r, 8, () as 



-r 2m+1 C 2n+1 cos 2m+2 0, 
1 + r 2m C 2n+1 cos 2m+1 6 sinO, 



C = ~^Zn 



r 2m+2 cos 2m+2 , 



2to+ 1 



r 2m C0f} 2 m ( 



With large enough r, the variables r and £ vary much 
slower than does. Therefore, we can replace cos 2m 9 
with its average c m defined by 



1 

2^ 



cos 2m Od0, 



and we obtain the following two equations: 

Jm+l a2ti+1 



r = -c„ 



2m+2 



2m+ 1 



(13) 
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The above two equations lead to 



_^dC=fr-^^|dr. 



P c m+ i r 



(14) 



With a function Z n (Q defined in 

dz n __ 2 C 2 " +1 

dC z n 

we have a conserved value H — (m + l)/3 _1 log_ffo + 
Z n determined by the initial condition. Here we have 
used the definition H n = r 2 /2 and the relation (2m + 
l)c m = 2(m + l)c m +i obtained from the integration by 
parts of Eq. (13). The function Z n (Q has the minimum 
value Z n (0) at ( = because dZ„/dC < if C < and 
dZ„/dC > if ( > 0. Therefore, we have the following 
inequality: 



Bu- 



rn 



1 



/3 



■logtfo <C 



(15) 



with a constant C . This inequality means that the energy 
of the system has the minimum and the maximum values 
(see Fig. 1), and that the system consequently loses its 
crgodicity. 

In order to confirm our arguments, we study three ther- 
mostats, i.e., g = p(, g = p( 3 , and g = p 3 (. All the 
thermostats are coupled with the one-dimensional har- 
monic oscillator and the inverse temperature [3 is set to 
be 1.0. For the relaxation time, we study three cases 
t = 5, 10, and 50. The initial condition is set to be 
(p,q,Q = (1.1,1-1,0). This condition gives C = 1.019 
for g — p( and g — p( 3 , and C — 0.829 for g = p 3 (. 
From the inequality (15), we can estimate the range of 
the energy as follows: 



0.815 <H < 1.211 (g 
1.210 <H < 3.076 (g 



pC,9 = pC 3 ), 
p 3 0- 



(16) 
(17) 



Time evolutions of the systems were numerically calcu- 
lated by the fourth-order Rungc-Kutta method with the 
time step 0.005 and those of the energies are shown in 
Fig. 2. The ranges of the energies agree well with our 
theoretical estimation (16) and (17) for larger values of r. 

We have studied the ergodicity of the thermostat fam- 
ily based on deterministic and time-reversible dynamics. 
We have obtained the conserved value for the harmonic- 
oscillator system coupled with the single-variable ther- 
mostats. This conserved value causes the energy to be 
bounded, and the system consequently loses its ergodic- 
ity. We performed numerical simulations and have con- 
firmed our theoretical arguments. The conserved value 
exists in the system with the single additional variable, 
since it is generally impossible to make a separable form 
as Eq. (14) for the system with two or more additional 
variables. Therefore, the number of the additional de- 
grees of freedom is essential for the ergodicity of the sys- 
tem [15]. While we have studied the harmonic-oscillator 
system, it is straightforward to apply our arguments 
for similar systems such as H = p 2 /2 + q 2k /2k (k — 
1,2,3, •••)• 

We have given the general expression for the one- 
variable thermostats with the pseudo friction term g = 
p2m+i£2n+i p or casc n = 0, the thermostat can be 
regarded as a method controlling higher moments of the 
kinetic energy [11]. On the other hand, there are no clear 
physical interpretations for general cases n/0. Addi- 
tionally, a more general form of pseudo Hamiltonian is 
available [12], while we have assumed the quadratic form 
of the additional variable as in Eq. (7). Therefore, it 
should be one of the further issues to clarify the physical 
meanings of general thermostats. 
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